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SECTION  I 


INTRODUCTION 

The  new  generation  of  supercomputers  (Cyber  205,  Cray  X-MP,  Cray  2,  and 
so  on)  heralds  In  a  new  era  of  research  for  aerodynamic  analysis  and  design. 

The  highly  complex  flow  field  about  a  weapon/store  configuration  on  current 
and  future  tactical  fighter  aircraft  can  now  be  modeled  and  examined  by  the 
computer.  Current  aerodynamic  research  at  the  Air  Force  Armament  Labora¬ 
tory's  Computational  Fluid  Dynamics  Section  is  aimed  at  solving  the  fluid 
flow  governing  equations  using  innovative  approximation  techniques.  These 
partial  differential  equations  form  a  nonlinear  system  of  equations  to  be 
solved  for  unknown  pressures,  densities,  temperatures,  and  velocities  that 
will,  in  turn,  yield  aerodynamic  characteristics  for  a  selected  weapon/store 
configuration  in  a  wide  range  of  flight  conditions. 

To  obtain  a  complete  understanding  of  the  flow  field  about  a  configura¬ 
tion,  we  must  solve  the  Navier-Stokes  equations.  One  of  the  first  steps  in 
developing  the  computer  codes  needed  to  solve  these  equations  is  to  do  an 
analysis  of  the  numerical  technique  selected  to  model  the  equations  of 
interest.  This  analysis  illustrates  the  procedures  used  to  evaluate  a 
numerical  technique's  ability  to  model  the  original  equations.  The  analysis 
yields  stability  criteria  and  limits  of  the  numerical  finite-difference 
technique. 

This  report  will  serve  as  an  instructional  aid  for  that  analysis.  A 
simple  model  equation  with  properties  characteristic  of  the  more  complex 
equations  of  fluid  flow  will  be  used  in  this  step-by-step  analysis.  The 
viscous  Burgers  equation  is  a  highly  useful  model  for  this  purpose  as 
explained  in  the  next  section.  The  method  of  MacCormack  will  be  the 
technique  employed  for  the  analysis. 

1.  BURGERS  EQUATION 

Within  the  field  of  computational  fluid  dynamics,  there  are  a  number  of 
simple  linear  problems  that  serve  as  vehicles  for  various  finite-difference 
methods.  Unfortunately,  the  usual  fluid  mechanics  problems  is  highly  nonlinear. 


A  single  equation  that  could  serve  as  a  nonlinear  analog  of  the  fluid 
mechanics  equations  would  be  very  useful.  This  single  equation  must  have  terms 
which  closely  duplicate  the  physical  properties  of  the  fluid  equations,  i.e., 
the  model  equation  should  have  a  convective  term,  a  diffusive  or  dissipative 
term,  and  a  time-dependent  term.  Burgers  (Reference  1)  introduced  a  simple 
nonlinear  equation  which  meets  these  requirements 
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Equation  1  is  a  parabolic  partial  differential  equation  which  can  serve  as  a 
model  equation  for  the  boundary-layer  equations,  the  parabolized  Navier- 
Stokes  equations,  and  the  complete  Navier-Stokes  equations. 


Burgers  equation  has  exact  analytical  solutions  for  certain  boundary 
and  initial  conditions.  These  exact  solutions  are  useful  when  comparing 
finite-difference  methods.  The  exact  steady  state  solution,  [i.e., 
lim  At  u(x,t)],  of  Equation  1  for  the  boundary  conditions 
t  — *  00 


u(O.t)  *  1  (2) 

u(L,t)  «  0  (3) 


is  given  by 


where 


(  1  +-  »  ip'  [  TT  R#i  (  Y  -  1  )  1  J 
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and  u  is  a  solution  of  the  equation 


For  simplicity,  the  linearized  Burgers  equation 
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(7) 


is  often  used  in  place  of  Equation  2.  In  this  analysis,  we  will  use  the 
nonlinear  case  for  the  computer  solution  and  the  linear  case  for  the 
analysis  in  Section  II. 


2.  MACCORMACK'S  METHOD 

The  MacCormack  method  (Reference  2)  for  the  complete  Burgers  Equation  1 
is 
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where 


F  *  0.5u2 


which  Is  second-order  accurate  In  both  time  and  space.  In  this  version  of 
the  MacCormack  scheme,  a  forward  difference  is  employed  in  the  predictor 
step  and  a  backward  difference  is  used  in  the  corrector  step  for  the 
convective  terms.  The  alternate  version  is  vice  versa  and  both  variants  are 
second-order  accurate.  However,  the  best  resolution  of  discontinuities 
occurs  when  the  difference  In  the  predictor  is  in  the  direction  of  propaga¬ 
tion  of  the  discontinuity  (Reference  3).  Since  we  will  be  giving  the 
discontinuities  a  positive  direction  as  initial  conditions,  [i.e.,  u(0,t)  *  1], 
we  will  use  the  version  indicated  above. 
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SECTION  II 
ANALYSIS 

In  this  section,  we  will  be  examining  the  features  of  the  linear  viscous 
Burgers  equation  (Equation  7).  The  MacCormack  method  for  Equation  7  is 
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The  first  step  will  be  to  combine  these  two  equations  into  one.  Let 
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To  provide  a  quick  check  of  this  tedious  combination,  note  that  the 
MacCormack  method  for  the  inviscid  Burgers  equation  is  equivalent  to  the 
Lax-Wendroff  scheme  (Reference  4).  Therefore,  if  we  set  «  =  0,  thereby 
removing  the  viscous  terms,  we  obtain 


2ur=  2u“i  -  *(«?♦■- 


The  next  step  In  the  analysis  will  be  to  check  consistency. 

1.  CONSISTENCY 

A  finite-difference  equation  (FDE)  is  said  to  be  consistent  with  the 
partial  differential  equation  (PDE)  It  models  if 

Llm  TE  *  0 
ax  -*•  0 
At  -*•  0 

where  TE  is  the  truncation  error  and  is  found  by  subtracting  the  original 
partial  differential  equation  from  the  finite  differential  equation. 

The  first  step  in  evaluating  the  consistency  of  the  finite  differential 
equations  obtained  is  to  express 

n+l  n  n  n  .  n 

UJ  •  Vl-  Vl-  “j+2'  “d  “j-2 

in  terms  of  Uj  an(j  ^ts  derivatives  by  using  Taylor  series  expansions,  i.e.. 


The  next  step  is  to  substitute  Equations  13  through  17  into  the  Equation  12. 
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Finally,  subtract  the  original  PDE  from  the  FDE.  Repeating  Equation  7  for 


reference  and  plugging  In  g  -  and  a 

Ax 


as  well  as  dividing 


through  by  At  yields: 
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Subtracting  Equation  7  from  Equation  19  yields  the  TE. 
that  Urn  TE  *  0  by  observing  that  the  ax  and  at  terms 


It  is  readily  apparent 


appear  only  In  the  numerators  of  all  terms.  Therefore,  the  FDE  given  by  the 
MacCormack's  method  for  the  linear  viscous  Burgers  equation  is  uncondition¬ 
ally  consistent. 


2.  STABILITY 

The  linear  viscous  Burgers  equation  will  also  be  used  for  this  portion 
of  the  analysis  of  the  MacCormack  FOE  so  that  the  Fourier  method  (or  Von 
Neumann  method)  of  stability  analysis  can  be  used. 

In  the  Fourier  or  Von  Neumann's  method,  the  numerical  stability  of  an 
FDE  Is  analyzed  by  introducing  a  disturbance  into  the  numerical  solution  at 
every  grid  point  In  the  spatial  domain  at  some  arbitrary  time  level.  That 
disturbance  is  then  expanded  into  a  Fourier  series  and  each  Fourier  compo¬ 
nent  of  that  series  is  analyzed  separately.  The  FOE  is  said  to  be  stable  if 
all  of  the  Fourier  components  do  not  grow  in  time  and  unstable  if  any  one  of 
the  Fourier  components  grows  In  time  (Reference  5). 

If  we  Introduce  a  disturbance  «  (which  represents  roundoff  errors)  into 
the  numerical  solution  at  every  grid  point  at  an  arbitrary  time  level  n 
using  Equation  12  we  get 
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Subtracting  Equation  12  from  Equation  20  yields 


h%+v2-  4K*i+*-i)  +  6?j 

Any  bounded  disturbance  «  on  IL  equally  spaced  grid  points  over  a  length  L 
can  be  represented  by  a  Fourier  series  made  up  of  a  finite  number  of  Fourier 
components:  B 
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where 
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Substituting  Equation  22  into  Equation  21  yields:  (Note:  also  changed  j 
subscripts  to  1) 
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Dropping  the  summation  sign  and  cancelling  the  common  exponential  factor 

Ik.x,  Ik, Ax 


,Ikjxi  yields: 
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Dividing  through  by  Aj  and  noting  that  a“+1/a“  -  Gj  (the  amplification 
factor)  yields 

(25) 
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-l-|/Sa[  2  (  cosK.dx  -I-  IsinKjAz)  -  2(cosKjdx  -  IsinK^x) 

-f-(cos2KjAx  -  I  s  i  n2KjA.x:)  -  (cos2Kjd.x’  +•  Isin2K.dx)] 
+^az[  (  cos2KjA.r  +-  Isin2K^x)  +■  (cos2^ax  -  I  s  i  n  2Kj ax  ) 

—  4-(cosK.  dx  -I-  IsinK.Ax)  -  4(cosKjd.x  -  IsinKdx)  +  6] 


or,  let  >*  kjax, 


Gj  *  1  ~2~P  (2Isinr) 

+  (*  Hfl[2(cos?  -  1)]  (26) 

•f j-]8o{ 41  s  i n 7  -  2Isin2?  ) 

+-|-a2(  2  c  o  s  2  7  ~  8cos7  +■  6) 

Using  the  following  trigonometric  identities 
cos2  y  *  2cos2y  -  1  and  sin2v  ■  2sinycosY 
we  obtain 

Gj  ~  l1  +  (2a  -»-02)  (cos  7  -  1)  -f  2  a2(  cos  7  -  l  )*  ]  ^ 

+■  I  |  -0  s  i  n.7[  2a  (cos7  -  1)  +■  l]j 

Now,  if  |Gj|^l,  the  disturbance  «  will  not  grow.  Therefore,  if 
lEquation  27|  <1 
or  (square  both  sides) 

lEquation  27|2<12 
Let  A  ■  cos  7-  1,  then 

[1  MW)A  **•«*•]• 

i"[  ~{i s  i  n.7  (  2aA  +-  1 )  ]*  <_  l 

Since  it  is  not  easy  to  obtain  simple  stability  criteria  with  the  three 
variables  in  Equation  28, use  of  a  simple  computer  program  can  determine 
these  criteria.  This  is  program  1  In  the  Program  Section.  The  first 
step  will  be  to  set  <*■  0  and  determine  the  limits  on  /?.  Remember  that 


Therefore,  when  «  *  0,  the  equation  reduces  to  the  inviscid  Burgers  equation. 
The  MacCormack  method  applied  to  the  Inviscid  Burgers  equation  Is  the  same 
as  the  Lax-Wendroff  method  mentioned  previously  and,  hence,  the  stability 
criteria  are  the  same,  l.e.,  q  £  1.  The  results  are  depicted  graphically 
in  Figure  1.  Note  that  when  / 8  »  1.1,  the  resulting  plot  fell  outside  the 
unit  circle. 

When  ft*  0,  the  equation  becomes  the  one-dimensional  heat  equation.  The 
closest  explicit  FDE  that  represents  this  case  Is  the  simple  explicit  finite 
differential  equation 


The  amplification  factor  for  this  scheme  is 

G  =  1  +  2 a(  cosr  -  1  ) 

while  for  FDE  (Equation  27) 

G  =  1  -f-  2«  ( co  s  7  -  1 )  +  2a^(cos7 

The  additional  term  does  not  affect  the  stability  criteria  since  both  are 
stable  when  «<  Figures  2  through  4  show  graphically  «  at  .25,  .5, 
and  .6.  Notice  that  at  «*  »  .6,  Instability  occurs. 

The  final  phase  of  the  stability  analysis  Is  to  try  to  determine  a 
relationship  between  <*  and  p  when  they  take  on  values  within  their  respec¬ 
tive  stability  ranges  (l.e.,  P<1  and«£y).  To  accomplish  this  task,  set  « 
at  a  fixed  value  and  vary  p  until  the  stability  limit  Is  reached.  Set  «  ».25 
and  .5  and  then  show  how  values  greater  than  .5  are  unstable  no  matter  what 
value  0  takes  on. 

First,  let  «*  .25  and  run  /Sat  0,  .25,  .75,  .8,  and  1  In  Figures  5(a) 
and  (b).  The  p  limit  is  Isolated  at  .86  in  Figures  6-9.  Next,  examine  «* 

.5  at  p  *  .5,  .75,  .9,  1,  and  2  In  Figures  10(a)  and  (b)»  The  ft  limit  is 
Isolated  at  .96  In  Figures  11  and  12.  Finally,  let  «  *  .75  and  run  p  at  .15 
and  .25  In  Figure  13.  Since  a  >  .5  In  this  case,  there  is  no  value  of  ft 
that  yields  stability. 


(30) 

-  .  i2  <31> 


The  final  results  are  tabulated  below.  These  values  will  be  used  as 
limits  in  the  rest  of  the  analysis. 


0.00 

0.25 

0.50 

0.50 


P 

1.00 

0.86 

0.92 

0.00 


-  Boundary  Limit 


-  Boundary  Limit 


An  exponential  curve  fit  for  the  values  of  «  and  p  between  their  limits 
(l.e.,  «*  .25  and  .5  while  .86  and  .92)  resulted  in  the  following 
relationship  for  stability 

p  £  0 . 7  7  0  4  e°  44a  (32) 

To  check  this,  plug  In  <*  ■  .3  and  .4  which  yields  p  *  0.8791  and  0.9187, 
respectively.  For  a  safety  margin,  round  down  to  two  significant  figures. 
Figures  14  and  15  are  plots  of  «  ■  .3  for  p  *  .88  and  .87.  Notice  in  Figure 
15  that  the  rounding  down  provided  an  appropriate  margin  of  safety  needed. 
Figures  16  and  17  are  plots  of  «  *  .4  for  p  *  .92  and  .91.  From  these  four 
plots,  one  can  see  that  Equation  32  can  be  used  confidently  for  the  relation¬ 
ship  between  «  and  p  within  their  respective  stability  ranges  not  to  Include 
the  boundaries. 

3.  CONVERGENCE 

For  this  section  of  the  analysis,  we  will  use  the  following  discussion 
from  Anderson  (Reference  3): 

"Generally,  we  find  that  a  consistent,  stable  scheme  is  convergent.  Con¬ 
vergence  here  means  that  the  solution  to  the  finite  differential  equation 
approaches  the  true  solution  to  the  partial  differential  equations  having  the 
same  Initial  and  boundary  conditions  as  the  mesh  is  refined.  A  proof  of  this 
is  available  for  initial  value  (marching)  problems  governed  by  linear  PDEs. 
The  theorem,  due  to  Lax,  is  stated  here  without  proof  (Reference  6):" 


"Lax's  Equivalence  Theorem.  Given  a  properly 
posed  Initial  value  problem  and  a  finite- 
difference  approximation  to  it  that  satisfies 
the  consistency  condition,  stability  is  the 
necessary  and  sufficient  condition  for  convergence." 


Since  we  have  shown  consistency  and,  if  we  stay  within  the  limits  of 
stability,  then  convergence  is  assured. 

4.  CONSERVATIVE  PROPERTY 

Consider  the  POE  with  the  following  boundary  conditions: 


iH.  +  c  — 

3t  3x 


32u 

3x2 


or 


JT  "K77(cu-yf3=0 


(7) 


(33) 


and  remember  c  and  are  constants.  The  BCs  are 
u(0,t)  -  Ux 
u(L,t)  *  0 

Equation  33  describes  the  following  conservation  principle  over  any 
control  volume  with  boundaries  at  x^  and  x^. 


Discretize  the  domain  as  shown  below 


i-  1  I  2  3  4  5  IL-1 1  IL 

1*1  1/2  and  at  1  *  IL  -  1/2  are  the  boundary  for  the  control  volume 
enclosing  all  interior  points 


Remember  that  MacCormack's  method  for  the  FDE  is  a  two-step  technique 
repeated  here 


Predictor: 


w+1 

u . 

j 


-  c-S-{u;M  -  up 


+  »-a 8ui  +  u?-.) 


(Ax) 


j+l 


(10) 


Corrector: 


u 


n*-l  _ 


-tK  + «; 


n+l  ,  n+t 


~  C 


(Hi 


A  t  t  ni-i  _  r 
+■  (u m  ~  2U, 


m-i 


n+  1 


(4i) 


1H 


‘jfl 


n+1  \  1 

Uj-1  )  J 


Considering  the  corrector  step  only  (without  plugging  in  the  predictor 
step's  algorithm) 

(35) 


2un+l  —  u“- 

J _ *  J 

At 


nTT  nTI 
Ui-1- 

c  — ZTx -  + 


nTT 

Vi 


j 


d+T 

Ui-1 


(Ax)- 


or 


«r-  u‘ 


At 


nTT 


pCU,.,  [I  .  »-»n+l 

=  1x5  “(Ax )a  J  J  ~ 


tHTI 


(36) 


rC_U, 

Ux 


for  j  «  2,  3,  _ IL-1 

Now,  sum  the  FDEs  given  by  Equation  36  over  the  control  volume  enclosing  the 
Interior  grid  points. 


'  .1  =  [CU1  “uT<U2'  “l)  ]  "  [  C  U8  -  U3  -  “a)] 


+  K  -iT^'h  - 


dTT 

»a>]  “ 


”  [cu.i.TT&K~  T 


nTT  ^TI 

=  tcui-^u2-  ui)J  ~  t  cu,,-,-^<un-  u, ,-,)  1 

Now,  remember  the  Predictor  step 

»r- »?-  (■*♦,- «?)  (1 

+  (ui«  -  2<  ^  cj 


Therefore,  LHS  of  Equation  37  becomes 


u"+1-  liT 


2  £  ^ i 

j-2  it 

with  the  RHS  of  the  Predictor  step  taken  to  the  RHS  of  this  equation. 
Dividing  through  by  2  yields: 


B+-1  D 
U:  -  Us 


■  -  [4*'  aif 

*  >]"  [‘4  vfi 


I l-l  (1-1 

L  (LHS)ix  =  2  { RHS ) lx 

i*2  j *2 


AMUUUUMQWML 


•/>.- '/.V 


Therefore, 


«  ,  i+i 

«»  +  U| 


»  ,  nt-1 

Vlii-i 


£(“  2+  C‘-  <-  O] 

-  2ta(“?,+  a",.,-  )] 


=  A  -  B 

Note  that  A  and  B  "sort  of"  represent  the  fluxes  at  i  =  1  1/2  and  at  i  =  I L 
-  1/2  (which  correspond  to  the  control  surfaces  of  the  control  volume).  It 
appears  that  the  convective  terms  of  A  and  B  are  questionable  whereas  the 
dissipative  terms  are  legitimate.  In  the  dissipative  terms 


for  example. 


we  obtain  an  averaging  of  the  predicted  and  corrected  values.  This  is  as 
expected  from  a  close  examination  of  this  two-step  method.  In  the  convective 
terms,  however,  note  the  forward  differencing  in  the  predictor  step  and 
backward  differencing  in  the  corrector  step.  For  this  reason,  within  the 
convective  terms,  uj  and  u^.j  are  from  the  predictor  step  while  and  uil 
are  from  the  corrector  step.  It  is  difficult  to  say  which  terms,  convective 
or  dissipative,  are  better  models  for  the  conservation  principle  in  this  FOE. 

The  averaging  within  the  dissipative  terms  may  smooth  things  out  too  much  or 
it  may  be  just  the  smoothing  needed.  The  mixing  of  steps  in  the  convective 
terms  may  leave  information  behind  or  add  data  too  soon  (similar  to  sinks  and 
sources).  In  conclusion,  MacCormack's  FDE  possesses  the  conservation  principle. 
As  with  all  else  in  numerical  analysis,  as  ax  and  at  approach  their  ideal  values 
within  stability  limits,  this  conservation  property  of  the  FDE  gets  stronger  and 
stronger. 

5.  TRANSPORTIVE  PROPERTY  AND  TRANSPORTIVE  ERROR 

For  fluid  flow,  heat  transfer,  and  combustion  problems,  a  disturbance 
at  any  location  in  the  flow  field  can  spread  to  other  locations  by  these 
mechanisms: 


a)  Convection  or  advectlon  (due  to  fluid  motion) 

b)  Diffusion  (due  to  random  molecular  motion) 

c)  Pressure  waves 


For  subsonic  flows,  diffusion  and  pressure  waves  can  spread  a  distur¬ 
bance  In  every  direction.  Convection,  however,  can  only  transport  the 
disturbance  in  one  direction;  the  direction  of  the  fluid  velocity. 

Since  the  linear  viscous  Burgers  equation  has  characteristics  similar 
to  the  viscous  incompressible  flow  equations,  we  expect  convection  and 
diffusion  transportive  properties,  the  MacCormack  method  FDE  will  possess 
the  transportive  property  if  it  has  these  properties. 

Consider  Equation  7,  the  linear  viscous  Burgers  equation 


(5 1 . 
H 


(7) 


We  have  seen  that  the  FDE  given  by  the  two-step  MacCormack  method  is 

3u’*'=  2u"  -  u$_,) 

+•  2“(“%i  -2u"  +  Vi1 
+  pz(u“+1-  Su1;  +  «1-,)  (12) 

+  2uf_,  +  u*„) 

+  4u"tl+  Ou,"-  4u“_,  +  «».,) 


where  a  -  y  -~r  and  B  ■  c  — 

(Ax)2  Ax 

Suppose  that  at  time  level  n,  u  is  zero  everywhere.  Now,  consider  a  distur¬ 
bance  y  >  0  at  grid  point  j  and  time  level  n  and  examine  how  this  disturbance 
spreads  to  the  neighboring  grid  points. 
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.u. 


2(un+1  -7)  =  (2a+p2)(-2  7)  4-  a2 (  6  7 ) 

at  j  (40) 

u°+1  =  7  (  1  +  3a2-  2a  -/?2) 


2<  vJ)  =  0(7)  +  (2a  +  P  c){y)  +  /2a  (-27)  +  aV*7) 

at  j+1 


at  j+2 


at  j-1 


at  j-2 


u;  +  i  -7(20  +2/2  -h  a  -f$ a-2  a  ) 
2<0  =  M?)  +  <#(7) 


VV  =  M 


V,'  =  ''(-fP+jf+a  +  pa  -  2  a2  ) 

V2  =  («2- Pa  ) 


,n+ 1  _  j_ 

~  2 

1  0  ,  1  a2 


(42) 


(43) 


(44) 


Recall  that,  for  numerical  stability,  we  found  these  two  criteria 

a  «  0,  g  <_  1, 

6-0,  a  <_  1/2 

For  a  *  0  at  time  level  n+1,  we  get 


j-2  j-1  j  j+i  j+2  x 


uj = 

7(1  -f) 

(45) 

Vl’ 

ly(P\P 

) 

(46) 

Uj+2  = 

0 

(47) 

Ui-1  = 

1  /  2 
f7  (P  -P 

) 

(48) 

Uj-2  = 

0 

(49) 

P  =  0,  Uj 

=  7 

— fc- 

1  _ 

V>  = 

Ui+2=  U,-l  = 

uJ_2=  0 
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(41) 


■  ui  =  T^  “,+1- fr .  v,= 


-v 


j-2  j-1  j  i*i  i+2  * 


f  ■  1  '  “j  -  V.=  ui-2=  V,-  0 


j-2  j-1  j  j+1  j+2  x 


Vl  =  1 


This  is  a  very  interesting  result  because  6  -  c  ^is  the  convection  term 


and  here  we  have  perfect  convection  at£  =  1. 


For  j3  ■  0  at  time  level  n+1,  we  get 


Uj  =  7  (3a  -  2a  +  1) 


«jfi*  r(«  -  2«2 ) 


Ui  +2  =  ) 


Uj_!=  7  (  a  -  2a2 ) 


vt-  ?(K) 


a  =  0  ,  n  i  =  7 


j-2  j-1  j  j+1  j+2  x 


Uifl=  Vl=  U  j  +2=  Uj-2=  0 


»v* 


j-2  j-1  j  j+1  j*2  x 


To  truly  represent  the  physics  In  the  process,  the  uj+2  and  uj_2  terms  must 
be  less  than  our  equal  to  the  uj+i  and  Uj_i  terms.  Therefore, 

1  2  2 

ja  <a-  za 

,  (55) 

fa<  1 

,  2  At  -  2 

all  or 


V 
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The  results  for  p  *  0  and  these  figures  tell  us 

1.  For  all  ««  [0,1/2],  the  disturbance  spreads  in  every  direction  (+ 
and  -  x-axis).  This  is  consistent  with  the  transport  characteristics  given 
by  Equation  2. 

2.  For  a  disturbance  of  the  form  i  (x-x.n  where  5(x-x.j)  is  the  delta 

,  At 

function, v  2  should  be  less  than  2/5  in  order  to  obtain  physically 
(Ax) 

reasonable  results  for  the  first  time  step. 

An  investigation  of  the  transportive  property  for  the  «  and  P  terms 
varying  within  their  respective  stability  ranges  is  difficult  analytically 
and  a  numerical  result  is  presented.  Please  note  that  the  program  used  to 


20 


generate  these  figures  models  the  original  nonlinear  viscous  equation  rather 
than  the  linear  equation.  We  give  the  domain  a  disturbance  at  u(15)  of  1  at 
time  level  0.  At  time  step  1  (Figure  18),  we  see  the  convective  property  in 
that  u(17)  is  0.1875,  and  we  see  the  diffusion  property  in  u(14),  u (16)  and 
in  u(17).  At  time  step  30  (Figure  19),  we  see  a  history  of  the  disturbances 
passing  through  the  discretized  domain.  Since  the  value  u  takes  on  between 
u(15)  and  u(29)  is  greater  than  Its  counterparts'  value  between  u(2)  and 
u(14)  at  this  time  step,  coupled  with  the  smoothing  out  of  the  disturbance, 
we  car.  say  that  the  FOE  has  the  same  transportive  properties  as  the  PDE. 
(Note:  u(l)  *  u(30)  *  0  are  fixed  boundary  conditions  within  the  code.) 

One  thing  to  point  out,  though,  is  that  we  have  been  looking  at  ;  =  .5  in  * 
Figures  18  and  19.  If  we  set  P-  .4,  the  limiting  value  for  a  true  repre¬ 
sentation  of  physics,  we  get  the  results  shown  in  Figures  20  and  21.  Notice 
how  these  plots  display  the  transportive  properties  in  a  more  realistic 
manner.  For  this  reason,  p  <  .4  is  recommended. 

6.  DISSIPATION  ERROR 

Diffusion  spreads  a  disturbance  in  every  direction  by  random  molecular 
motion.  As  diffusion  spreads  a  disturbance  in  every  direction,  the  distur¬ 
bance  is  smoothed  out  over  an  increasingly  larger  region  of  space,  thus 
reducing  spatial  gradients  and  lowering  the  magnitude  of  the  disturbance. 
This  effect  is  called  dissipation. 

Diffusion  is  mathematically  described  by  even-order  derivatives,  i.e., 
c  92nu  .  c  32nu  n*l,2,...  Specific  examples  include  yAV  in  the 

n  3x^2n  n  3x*3x” 

momentum  equation  which  represents  the  diffusion  of  momentum  and  kv2y  -jn  the 
thermal  energy  equation  which  represents  the  diffusion  of  thermal  energy. 

If  all  of  the  coefficients  of  the  even-order  spatial  derivatives  in  a  PDE 
are  zero,  then  that  POE  does  not  have  any  dissipative  characteristics. 

In  order  for  an  FDE  to  have  the  same  dissipative  characteristics  as 
the  PDE  that  it  is  to  represent,  the  coefficients  of  the  even-order  spatial 
derivatives  In  the  modified  equation  (for  the  FDE)  must  be  identical  to 
the  corresponding  coefficients  in  the  POE.  The  modified  equation  is  the 


POE  which  is  actually  solved  when  a  finite-difference  method  is  applied  to 
a  POE.  If  the  corresponding  coefficients  are  different,  then  the  solutions 
of  the  FDE  will  contain  dissipation  errors. 


The  dissipation  errors  in  the  solutions  of  FDEs  are  usually  examined  by 
the  Von  Neumann's  method.  The  dissipation  error  for  each  Fourier  component 
of  a  disturbance  after  n  time  steps  is  given  by 

cl gPDeI (n)  -  |G](n)]A°  (56) 

where 


Gpde=  amplification  factor  of  the  PDE 
G  ■  amplification  factor  of  the  FDE 

A0  *  initial  amplitude  of  the  Fourier  component 

|*J  *  modulus  of  * 

Dissipation  error  is  given  by  Equation  56  because  the  modulus  of  the  ampli¬ 
fication  factor  depends  solely  on  the  even-order  spatial  derivatives  when 
the  highest  time-derivative  is  first  order. 


Once  again,  consider  the  linear  viscous  Burgers  Equation 

77  +•  c  ~  *  n  c  and  m  are  real  constants  (7) 

6t  <J*  ^  d%7 

This  equation  obviously  has  dissipative  characteristics  since  the  coeffi¬ 
cient  of  the  even-order  spatial  derivative  is  not  zero.  To  find  the  ampli¬ 
fication  factor  of  Equation  2,  substitute  a  Fourier  component 


uj  .  e 

where 

kj  *  wave  number 

I  *  V^T 


(57) 


«  *  unknown  complex  number  (a+bI)[Equation  7  will  be  used  to  determine 

a,  l.e.,  only  disturbances  that  are  also  solutions  of  Equation  7  are 
considered] 


of  the  disturbance  u  ~  ^  ui 

j  1 


into  Equation  7  and  obtain 


c  I  k^**1  e1*1  =  M(lk/  e  °'1  e1  ^ 

Since  I2  *  -1  and  eliminating  uj  ■  )  tg  I  k j  x  t  we  ge^ 

Oj  -  -iik‘-  elk,  ,  i.e.  ,  a  =  -/ik‘  and  b  =  -ck,  (58) 
By  using  Equation  58,  the  amplification  factor  for  Equation  7  is 


GPDE  =  eIfrPDE 


Since 


GPDE  =Igpde»  e  PDE 


Ir  I  =  At 

I^Dnc*  e 


and,  therefore, 


|GpnCl  S  1 


We  see  the  disturbance  is  dissipated  as  expected.  Let  us  now  evaluate 
the  dissipation  error  for  our  FDE  repeated  below 

2«r- 2u"  -  *(«•♦!-  “I-.) 

+  2a  ( u“f,  -2u“  +  (12) 

+  da(2Uj*,-  2“°-i  +  u°-2~  “%s) 

+  a*(u"jt2-  4u“>1+  8u?  -  4u“_j  +  u‘.2  ) 


*1  ,»*  -'» 


'll 


J 


The  modified  equation  for  Equation  12  is 

+  c  |E  _  -  -  f  (Ax2-c2At2) 

3t  3X  3x2  6  3*3 


cV  2  2a  2,  A  v  ,  2a  2  aA  A 


c  At  /  2  2.  2.  ,  V  /  2.  /  ax  s  d  U 
-  -=^-  (Ax  -c  At  )  +  -J  (c  At - g“)  — 2  •  *  * 

Since  all  the  coefficients  of  the  even-order  spatial  derivatives  in  the 
modified  equation  do  not  match  all  the  corresponding  coefficients  in 
Equation  7,  there  is  a  dissipation  error  in  the  solution.  The  modified 
equation  is  said  to  have  fourth-order  dissipation  because  the  second-order 
dissipative  term  coefficient  matches  the  dissipative  term  coefficient  of 
the  PDE. 

The  amplification  factor  for  Equation  7  is  given  be  the  Von  Neumann's 
method  (from  the  stability  analysis  done  earlier)  as 

0^  —  (l  +  (2  a  -f-)S2)  (cos  7  -  1)  +■  2a2(cos  7  -  1)  ]  (27) 

+  Ij-0sin7[2a(cos7  -  1)  +  1  ]  | 


where 


y  -  yx 


If  A  *  (cos  r-  1),  then  the  modulus  of  the  amplification  factor  is 
|G|  =  |  [  1  +  (2a  -f- Z?2 )  A  +-  2a2  A*] 

i 

+  [  -0  s  i  n  7  (  2aA  +  1  )  ff 


The  dissipation  error  for  each  Fourier  component  of  the  solution  of 
Equation  7  after  n  time  steps  is  (assuming  initial  amplitude  of  unity) 

|«Posl,B>  -  H(n) 


Recall  that  iGpggl  *  e^j  At  and  >  *  kjAX*  Thus,  the 
following  can  be  written 

IG. 


s  AA. 
-T 


PDE  1 


JPDE 


-7  « 


Programs  2  and  3  in  the  program  section  generate  this  curve  and  the  curve 
for  | G |  in  the  figures  to  be  discussed  next. 


When  a  *  0,  the  result  is 
iGpDEl*  1 

and  so  the  dissipation  error  Is 
1  -  |G|(°) 

Note  that  when  a  =  0,  the  Lax-Wendroff  scheme  is  obtained.  Figure  22  is  the 
plot  for  the  dissipation  error  in  this  case  which  agrees  with  previous  work 
on  the  Lax-Wendroff  scheme  (Reference  4).  The  dissipation  error  is  small 
when  p  (also  known  as  the  Courant  number)  is  small,  not  just  when  r  is 
small.  This  indicates  that  when  p  is  small  (and  «  *  0),  the  dissipation  for 
every  Fourier  component  is  small. 


For  a  at  other  values  within  its  stability  limits,  the  following  must  be 
considered 

(lW(n)  - 

where  iGpQglls  defined  above.  First,  set  «*  .25  and  run  P  at  .25,  .5,  .75, 
and  .86  -  values  within  the  stability  limits  for  this  «.  Also  run  p  at  .9 
for  curiosity's  sake.  In  Figure  23,  the  dissipation  error  grows  as  r  in¬ 
creases  at  the  same  rate  for  all  values  of  P  examined.  For  larger  values  of 
j8,  the  error  grows  faster.  Next,  set  «  *  .5  and  run  P  at  .25,  .5,  .75,  .96, 
and  1  In  Figure  24.  The  results  are  roughly  the  same  as  at  <*  *  .25  in 
Figure  23. 


For  many  problems,  the  solution  depends  primarily  on  the  first  few  terms 
of  the  Fourier  series  which  correspond  to  Fourier  components  with  low  wave 
numbers  (i.e.,  long  wavelengths).  In  these  cases,  it  would  be  acceptable  to 


have  a  large  amount  of  dissipation  error  in  Fourier  components  with  large 
wave  numbers  (1.e.t  short  wavelengths)  since  these  terms  are  very  small  and 
do  not  contribute  significantly  to  the  solution. 

Furthermore,  since  Fourier  components  with  the  largest  wave  numbers 
are  the  least  stable,  it  is  often  desirable  to  have  more  dissipation  at 
the  larger  wave  numbers  to  assure  stability. 

FDEs  with 

|g|  <  1 

for  every  Fourier  component  are  said  to  be  dissipative  FDEs.  Note  that 
|G|^1  also  indicates  numerical  stability.  The  total  error  in  the  solutions 
of  FDEs  is  made  up,  in  general,  of  two  parts:  discretization  error  and 
stability  error.  Stability  error  is  small  for  stable  FDEs  since,  by 
definition,  disturbances  and  errors  cannot  grow  for  stable  FDEs.  Thus, 
discretization  error  accounts  for  most  of  the  total  error.  The  discretiza¬ 
tion  error  can  be  broken  up  into  two  parts:  the  dissipation  error  and  the 
dispersion  error.  Recall  that 

u"+l  -  Gun 

so  that  if  G  and  u  are  correct,  then  u  will  be  correct  and 

G  =  |G|  e1* 

The  error  in  |G|  gives  rise  to  dissipation  error  and  the  error  in  ♦  (the 
phase  angle)  gives  rise  to  dispersion  error.  Dispersion  error  will  be 
discussed  In  the  next  section. 

Since  | G I  Is  related  to  even-order  spatial  derivatives  and  *  is  related 
to  odd-order  spatial  derivatives,  dissipation  error  dominates  over  disper¬ 
sion  error  If  the  lowest-order  term  in  the  truncation  error  is  an  even-order 
spatial  derivative.  The  reverse  is  true  If  the  lowest-order  term  in  the 
truncation  error  is  an  odd-order  spatial  derivative.  For  the  MacCormack 
method  FOE  examined,  note  from  the  modified  equation  that  the  lowest-order 
spatial  derivative  is  of  odd-order.  Therefore,  dispersion  error  dominates. 


In  fluid  flow,  heat  transfer,  and  combustion  problems,  dissipation  error 
arises  from  differencing  the  convection  terms  (first-order  spatial  deriva¬ 
tives).  Dissipation  error  should  be  minimized  so  that  numerical  dissipation 
does  not  swamp  physical  dissipation.  From  the  coefficient  for  the  fourth- 
order  spatial  derivative  in  the  modified  equation  for  the  MacCormack  FDE,  a 
functional  relationship  can  be  derived  between  a  and  p  that  will  eliminate 
this  error  term  as  follows 


^(Ax2  -  c2At?)  =  -£(¥*-  c2AC) 


1  [f-fi 

a  = 

A  plot  of  this  relationship  is  shown  in  Figure  25. 


7.  PHASE  AND  DISPERSION  ERROR 

A  medium  Is  said  to  be  dispersive  if  the  phase  velocities  of  waves 
propagating  in  that  medium  depend  on  the  wavelengths  of  the  waves.  In  a 
dispersive  medium,  waves  with  different  wavelengths  will  separate  more  and 
more  from  each  other  In  space.  The  phase  velocity  of  a  sinusoidal  wave 

eat  eIkix  =  e*1  gMkjx  +  bi)  »  e*1  eIki*x+  I}*)  is 
Dispersion  is  mathematically  described  by  odd-order  spatial  derivatives, 

3  u 

1.e.,cn  2n+I  »  n«l,2,...  if  the  highest-order  time  derivative  is  first- 

order.  In  order  for  an  FDE  to  have  the  same  dispersive  characteristics  as 
the  PDE  that  It  is  to  represent,  the  coefficients  of  the  odd-order  spatial 
derivatives  in  the  modified  equation  (for  that  FDE)  must  be  identical  to  the 
corresponding  coefficients  In  the  PDE.  If  the  corresponding  coefficients 
are  different,  then  the  solutions  of  the  FDEs  would  contain  dispersive 
errors. 

Dispersion  error,  like  dissipation  error,  is  also  examined  by  the  Von 
Neumann's  method.  The  dispersion  error  for  each  Fourier  component  of  a 
disturbance  after  n  time  steps  is  given  by 


where 


MW*) 


(60) 


*PDE  a  phase  angle  of  the  amplification  factor  of  the  PDE 
*  *  phase  angle  of  the  amplification  factor  of  the  FDE 

The  relative  phase  shift  error  for  a  given  Fourier  component  after  a  one 
time  step  is  defined  as 

£ 

*m.  (61) 

If  the  relative  phase  shift  error  >1  for  a  given  Fourier  component,  then 
the  numerical  solution  for  that  Fourier  component  gives  a  wave  speed  >  the 
wave  speed  given  by  the  exact  solution.  This  is  called  a  leading  phase 
error.  If  the  relative  phase  shift  error  <  1,  then  the  numerical  solution 
gives  a  wave  speed  <  the  wave  speed  of  the  exact  solution.  This  is  called 
lagging  phase  error. 

Dispersion  error  is  given  by  Equations  60  and  61  because  the  phase  angl 
of  the  amplification  factor  depends  solely  on  the  odd-order  spatial  deriva¬ 
tives  when  the  highest  time  derivative  is  first-order. 


Once  again,  consider  the  linear  viscous  Burgers  equation 


h.  +  c  is.  c  n 

<5*  ^  6*: 


it 


(7) 


where  c  and  m  are  real  constants.  The  amplification  factor  for  Equation  2 
was  derived  in  the  previous  section  where  the  phase  angle  *pqe  is  * 


-ckj  At 


-  By,  B 


it 

ca5T  •  r 


kj  Ax 


The  phase  velocity  is 


c  and  is  the  same  for  every 


Fourier  component.  Thus,  Equation  7  does  not  possess  any  dispersive 
characteristics. 


Evaluation  of  the  dispersion  error  of  the  MacCormack  method  applied  to 


Equation  7  yields: 

n+l 


n+i  n  A  t  /  „  „  . 

uj  “  -  up 


Predictor: 


+  +  V.) 


Corrector:  „n+i  _  J_  r  „n 


ur1=^[ui  +  u“+1  -  cii(^+i-  n£) 


(10) 


(11) 


n  +  l 


2  v  “j+1  *uj«  T  )  ] 


„  At  i  n+l  i 

*  ^~2(ui+i  -  2uIit  + 

(Ax) 

The  modified  equation  for  the  FDE  is 

£  *ci£  -  „  il=  -  -  |(Ar‘  -  C*AtI)^~ 

Si  St  *  St*  6  *6  iJ 

-[^<Al‘-c‘At)  +  f  (c‘Af-  .  .  . 

Since  the  coefficients  of  the  odd-order  spatial  derivatives  in  the  modi¬ 
fied  equation  do  not  match  the  corresponding  coefficients  in  Equation  7, 
the  solutions  of  the  FDE  contain  dispersion  error  of  third-order.  The 
order  of  dispersion  Is  equal  to  the  order  of  the  lowest  order  odd-order 
spatial  derivative  with  nonzero  coefficients  excluding  the  first-order 
spatial  derivative.  In  general,  the  higher  the  order  of  dispersion,  the 
lower  the  dispersion  error. 


The  phase  angle  of  the  amplification  factor  for  Equation  7  is 


$  —  t  arr'  lm(G) 

Un  STTi) 


sin7( 2 a  A  +  1 ) 

FDE  1  +  (2a+^)A  +  2az# 


$  = 


(62) 


At  At  , 

a  ■  u  - r;  6  ■  c  — »  A  “  cosy-l ,y“k . Ax 

(Ax)2  1 


where 


The  dispersion  error  for  each  Fourier  component  of  the  solutions  of  the 
FDE  after  n  time  steps  is 

t  - 1  -f}  s  inf  (  2a  k  +  1 )  "1 

*ai1  l  +  (2a+^)A  +  2aZkZ  -•  (64) 

The  relative  phase  shift  error  for  each  Fourier  component  after  one  time 
step  is 


n(W*)  =  n[“  “ 


±_  = 
^PDE 


-fi  s in7 ( 2a  k  +  1  ) 
t  aif1  1  +  (2a+^)A  +  Sa2!? 

-  py 


(65) 


For  «  =  0,  or  the  Lax-Wendroff  scheme  for  the  linear  inviscid  Burgers 
equation,  a  plot  for  p  =  .25,  .5,  .75,  and  1  is  in  Figure  26.  The  solution 
here  has  predominantly  lagging  phase  error  except  for  large  wave  numbers  with  /?« 
(vC57  1.)*  For  0*1*  the  relative  phase  shift  error  increases  as  r  increases. 
From  Figure  22,  notice  that  when  «  *  0,  the  FDE  is  not  very  dissipative  for 
any  p  when  r  is  low  (Figure  22).  As  a  result,  dispersion  error  is  very 
important  when  <*  *  0.  For  «  =  .25,  a  plot  of  P  *  .25,  .5,  .75,  .86,  and  .9 
is  in  Figure  27.  Observe  that  the  solution  here  has  lagging  phase  error  for 
all  wave  numbers  when  Pt  (0,.5)  which  increases  as  r  increases.  For  pe 
(.5,.9),  small  lagging  phase  error  is  observed  for  y  up  to  60°.  When  y  >  60  , 
the  leading  phase  shift  error  takes  over  which  increases  as  y  increases.  At 
approximately  y>  110° ,  the  leading  error  decreases  as  T  increases.  For<*=  .5, 
a  plot  of  P=  .25,  .5,  .75,  .96,  and  1  is  in  Figure  28.  For  p*  (0,.5), 
lagging  phase  error  occurs.  For  Pe  (.5,  1),  leading  phase  error  occurs  which 
Increases  dramatically  as  y  increases.  The  same  conclusions  that  were  drawn 
from  Figure  22  for  the  «  *  0  case  can  be  drawn  from  Figures  23  and  24  for  these 
two  cases  since  the  general  trend  seems  to  be  the  same.  The  program  that 
generated  the  data  for  this  figure  is  in  the  Program  Listing  Section  (Program 
4). 


30 


SECTION  III 


COMPUTER  SOLUTION 


For  this  section  of  the  analysis 
will  be  solved 


6u> 

Ti 


+  Uj 


,  the  nonlinear  viscous  Burgers  equation 


(1) 


Equation  1  will  be  solved  numerically  by  rewriting  it  in  conservative  law 
form 


6  u 
6  t 


(66) 


The  conservative  law  form  of  PDEs  occurs  when  the  coefficients  of  the 
derivative  terms  are  either  constant,  or,  if  variable,  their  derivatives 
appear  nowhere  in  the  equation.  This  form  of  the  PDE  circumvents  some 
numerical  problems  where  discontinuities,  such  as  shock  waves,  show  up. 


The  MacCormack  method  for  this  PDE  is 


Predictor: 


n>l  n  Air  /!_  /!_  2vni 

Uj  =  Uj  -  sl  V2  u  )jM  “  (2  U  ).J 


(67) 


Corrector: 
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(68) 
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A  listing  of  the  exact  algorithm  used  to  implement  MacCormack's  method  is  in 
Program  5.  Within  the  program  were  options  for  boundary  condition  values 
which  remain  fixed  during  program  execution.  One  could  also  change  the 
starting  solution  as  was  demonstrated  In  the  Transportive  Property  section 
when  a  point  disturbance  was  Introduced  in  the  middle  of  the  discretized 
domain. 


simplified  the  analysis  of  other  parameters.  The  MR  parameter  in  the 
output  figures  represent  the  mesh  Reynolds  number  in  a  sense.  It  actually 
is  —  since  the  velocity  type  term  is  variable  here  rather  than  a  constant 

M  At  1  A, 

c.  This  MR  parameter  also  tells  about  P  =  v  2  which  is  just  (since-— 

(.Ax)  r®  Ax 

■  1).  The  domain  is  discretized  into  30  sections  and  the  time  parameter  can 
run  as  long  as  one  wishes.  The  parameter  m  is  the  only  variable  that  will 
be  used  here  besides  time.  Also,  u ( 1 )  =  1  and  u(30)  =  0  for  all  cases.  The 
stability  limits  and  other  results  from  Section  II  can  be  extrapolated  if  a 
small  safety  margin  is  used  (to  account  for  the  crossing  into  the  nonlinear 
domain  from  the  linear  realm). 

Begin  with  m  =  .1  which  yields  MR  *  1  and  hence  P  =  1.  The  /?  limit 
derived  earlier  for  the  linear  case  was  P  <  .5  so  instability  is  expected 
in  this  case.  Figures  29  and  30  show  the  first  two  time  steps  at  these 
values,  and  the  solution  goes  unstable  very  fast.  When  m  =  .4,  MR  =  .25, 
and  hence  p  =  4,  the  instability  grows  even  faster  as  can  be  seen  in  Figure 
31.  As  the  value  of  MR  approaches  2,  i.e.,  p  approaches  .5,  the  effects  of 
dispersion  error  lessens  and  a  more  stable  solution  evolves.  Figures  32 
through  37  show  the  histories  of  the  solution  as  p  takes  on  values  of  .67, 
.57,  .556,  .54,  .53,  and  .5.  Notice  in  Figure  37  also  that  the  solution  is 
very  smooth  and  resembles  the  exact  solution  with  dissipative  characteristics 
extremely  well. 

The  analysis  for  the  transportive  property  showed  that  for  the  linear 
FDE,  a  value  of  P  <  .4  models  the  physics  of  the  diffusion  process  in  a 
more  realistic  sense.  When  p  is  run  at  values  of  .4,  .38,  .36,  and  .33  in 
Figures  39  through  42,  it  can  be  seen  that  as  the  value  of  p  decreases  the 
amount  of  time  it  takes  for  the  dispersion  error  to  grow  decreases.  From 
this  small  test,  it  would  seem  that  keeping  P  between  .4  and  .5  would  yield 
the  most  robust  solutions. 


SECTION  IV 


CONCLUSIONS 

This  report  has  illustrated  the  procedures  used  to  evaluate  a  numerical 
technique's  ability  to  model  a  partial  differential  equation.  The 
MacCormack  numerical  technique,  as  applied  to  the  linear  viscous  Burgers 
equation,  was  used  for  the  analysis.  The  finite  difference  equation  from 
the  MacCormack  method  was  examined  for  consistency,  stability,  convergence, 
dissipation,  phase  and  dispersion.  A  computer  solution  using  the  MacCormack 
method  was  compared  to  the  analytical  results.  Valuable  stability  limits 
and  parameters  for  maintaining  accurate  physics  modeling  determined  in  the 
analytic  study  were  demonstrated  with  the  computer  solutions. 

The  procedures  documented  herein  are  applicable  to  a  wide  range  of 
problems  that  require  an  analytical  evaluation  of  a  numerical  technique 
before  actual  computer  implementation.  This  analytical  evaluation  is 
necessary  if  computer  results  are  to  be  used  with  confidence. 
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Figure  3.  Stability  for  heat  equation  at 
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Figure  4.  Instability  of  heat  equation  at  a> 
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Figure  13.  Instability  at  a  -  .75  and  B  -  .25  and  .15. 
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Figure  17.  Closeup  of 


BURGERS  EQUATION  SOLUTION  time  step 
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Figure  18.  Transport ive  property  for 


BURGERS  EQUATION  SOLUTION  time  step 
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Alpha  vs.  Beta  Relationship 
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Figure  33.  Instability  at  B  *  0.57 


CD  ®  Lf>  ®  T-  ®  m 

O  6}  O)  <S  O  00  S 
©  (6  ®  ®  O)  O 

6D  S  fi  ®  ®  S  ffi  O 
(11850)0(0  0)0 
O  O  O)  o  o  o)  O 
•  •  O)  •  ■  0)  • 


o^mo-s-com® 

LT)®tv-®®'*— ■«— oj 
O)CSO)03SO)SN 
0)00)000)00) 
O)  ®  0)  0  0  0)0  O) 
0)00)000)09) 
O)  •  O)  •  •  0)  •  O) 


U)r  nOWO)r 

t  cij  oj  n  oj  t  oi 

00  O  O  0-  TT  U)  O) 

O)  O  O  O  O  x-  cu 

0)000000 
0)000000 
O) . 


«3MJ)(*)MDNO 

no)r-or  smis 
oo'r-oorur-coroo 
ojcyor-o)0)t^-o 
ooo)0^-m'r-cs-o 

0)00)0)U)i—  0)0 
0)0)00C^-0)aiiVO 


o  o  o 


o  o  o 


CD  (2  Q  0  0  (S3  0  (2 


CO  n  n  n  h  n  n  a  n  m  n  n  n  n  n  n  n  ii  ti  ii  ii  n  n  it  ii  ii  ii  n  ii  ii  it 


©oo®cDift*-r'-co*-cneo-r-<x>ojcocncuO'r-*-oju>ooLfti-,rt'-*-oooD 

SOIOOOtOansanOOnOOKOlA^OiatOOnratNOOO 

m  m  m  m  m  ®  ©  a  ©  m  m  ©  o>  o>  m  m  a>  a>  ©  ©  ©  ru  ld  m  t'-  *-  r-  oo  © 

®®®SBOiK@ffi®@e>®0)©aRC!!0)5)iDffl!Ss»-r(BWnt'*® 
(OSWSIBftSBftOSSWOJSffiC^CDttOsOSISQKISOOWWOnS 
O  K  9  O  CS  S)  O)  S)  v:  C;i  S)  t?  Ol=  (I  O)  ®  W  O  O  5)  O  S)  CD  O  W  Lf)  00  O 

•  o)  •  •  o>  •  o>  o>  *o>o)  •  o>  o>  o>  o)  o>  o> . •  o>  o>  oo  ©  © 


•  b>  o>  • 


O)  O)  oo  ©  © 


CO  O  (Q  CO  ©  ^ 


©  CO  S)  O  (S 


nr 


HMiiMHiinHHnitnniinnnnniiiinnnniinnn 

<uco^Lft©t^ooo>®<»-(um'j-©©r'-ooa>®'«--ajcnTr©©r'-coc&© 

^->r-^T-^-^-T-T-»-'r-rufuaiojfurucuajaiairo 


Z3  Z3  3  =>  3  ©  Z3 


Ifto- 
£  «X 

-c~a 

-O  «  \ 


il  1 00C 
■>  i.Ui 

q  a 
£  « 

(.  30 

o  or 
(J  oo 

o  *><j 


£3W 


mmmmw 


Figure  35.  Stability  at  6  ■  0.54 
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Figure  36.  Stability  at  “  *  0.53 
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